/*************************************************************************
  > File Name: propagation.c
  > Author: wusihai
  > Mail: wusihai18@gmail.com 
  > Created Time: 三  1/29 12:50:10 2020
 ************************************************************************/
#include "rtmlib.h"
#include "slavertm.h"
#include <athread.h>
#define PI 3.14159265354

extern SLAVE_FUN(Supdate_Vel)(rtm_slave_parameter* p);
extern SLAVE_FUN(Supdate_T)(rtm_slave_parameter* p);
//#define SUNWAY
/**************************
 * 正向传播相关函数*
 * ***********************/
//加载震源
void add_souce(int sx,int sy,int sz,int it){
  if(it<wave_len){
	Txx[n3d_ext(sx,sy,sz)]+=wavelet[it];
	Tyy[n3d_ext(sx,sy,sz)]+=wavelet[it];
	Tzz[n3d_ext(sx,sy,sz)]+=wavelet[it];
  }
}
//更新速度场
void update_Vel(){
#ifdef SUNWAY
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;
  //合并应力分量
  for(i=0;i<SubNx;i++){
	for ( j = 0; j < SubNy; j++){
	  for ( k = 0; k < SubNz; k++){
		T[n3d_ext(i,j,k)]=Txx[n3d_ext(i,j,k)] +Tyy[n3d_ext(i,j,k)]+Tzz[n3d_ext(i,j,k)];
	  }
	}
  }

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float rho=1.0f;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;

  //参数输入
  rtm_slave_parameter parameter;

  parameter.PT=T;
  parameter.PTxx=Txx;
  parameter.PTyy=Tyy;
  parameter.PTzz=Tzz;
  parameter.PVx=Vx;
  parameter.PVy=Vy;
  parameter.PVz=Vz;
  parameter.Pvel=vp_sub;

  parameter.SubNx=SubNx;
  parameter.SubNy=SubNy;
  parameter.SubNz=SubNz;
  parameter.stX=stX;
  parameter.endX=endX;
  parameter.stY=stY;
  parameter.endY=endY;
  parameter.stZ=stZ;
  parameter.endZ=endZ;

  parameter.dt=dt;
  parameter.dh=dh;
  parameter.rho=rho;
  parameter.Vpmax=Vpmax;


  parameter.DivNx=DivNx;
  parameter.DivNy=DivNy;
  parameter.DivNz=DivNz;
  parameter.x_rank=x_rank;
  parameter.y_rank=y_rank;
  parameter.z_rank=z_rank;

  parameter.PMLX=PMLX;
  parameter.PMLY=PMLY;
  parameter.PMLZ=PMLZ;

  parameter.mgridnumx=mgridnumx;
  parameter.mgridnumy=mgridnumy;
  parameter.mgridnumz=mgridnumz;

  //启动从核
  athread_spawn(Supdate_Vel,&parameter);
  //回收从核
  athread_join();

#else
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float rho=1.0f;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;

  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for (l=0;l < 5;l++){
		  temp1 += StaCoef[l] * 
			( Txx[n3d_ext((i+1+l),j,k)] - Txx[n3d_ext((i-l),j,k)]
			  +Tyy[n3d_ext((i+1+l),j,k)] - Tyy[n3d_ext((i-l),j,k)]
			  +Tzz[n3d_ext((i+1+l),j,k)] - Tzz[n3d_ext((i-l),j,k)]);

		  temp2 += StaCoef[l] * 
			(Txx[n3d_ext(i,(j+1+l),k)] - Txx[n3d_ext(i,(j-l),k)]
			 +Tyy[n3d_ext(i,(j+1+l),k)] - Tyy[n3d_ext(i,(j-l),k)]
			 +Tzz[n3d_ext(i,(j+1+l),k)] - Tzz[n3d_ext(i,(j-l),k)]);

		  temp3 += StaCoef[l] * 
			(Txx[n3d_ext(i,j,(k+1+l))] - Txx[n3d_ext(i,j,(k-l))]
			 +Tyy[n3d_ext(i,j,(k+1+l))] - Tyy[n3d_ext(i,j,(k-l))]
			 +Tzz[n3d_ext(i,j,(k+1+l))] - Tzz[n3d_ext(i,j,(k-l))]);
		}
		Vx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*Vx[n3d_ext(i,j,k)] + dt / rho* (temp1) / dh);
		Vy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*Vy[n3d_ext(i,j,k)] + dt / rho* (temp2) / dh);
		Vz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*Vz[n3d_ext(i,j,k)] + dt / rho * (temp3) / dh);
	  }
	}
  }
#endif
}
//更新应力场
void update_T(){

#ifdef SUNWAY
 //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  //参数输入
  rtm_slave_parameter parameter;

  parameter.PT=T;
  parameter.PTxx=Txx;
  parameter.PTyy=Tyy;
  parameter.PTzz=Tzz;
  parameter.PVx=Vx;
  parameter.PVy=Vy;
  parameter.PVz=Vz;
  parameter.Pvel=vp_sub;

  parameter.SubNx=SubNx;
  parameter.SubNy=SubNy;
  parameter.SubNz=SubNz;
  parameter.stX=stX;
  parameter.endX=endX;
  parameter.stY=stY;
  parameter.endY=endY;
  parameter.stZ=stZ;
  parameter.endZ=endZ;

  parameter.dt=dt;
  parameter.dh=dh;
  //parameter.rho=rho;
  parameter.Vpmax=Vpmax;

  parameter.DivNx=DivNx;
  parameter.DivNy=DivNy;
  parameter.DivNz=DivNz;
  parameter.x_rank=x_rank;
  parameter.y_rank=y_rank;
  parameter.z_rank=z_rank;

  parameter.PMLX=PMLX;
  parameter.PMLY=PMLY;
  parameter.PMLZ=PMLZ;

  parameter.mgridnumx=mgridnumx;
  parameter.mgridnumy=mgridnumy;
  parameter.mgridnumz=mgridnumz;

  //启动从核
  athread_spawn(Supdate_T,&parameter);
  //回收从核
  athread_join();

#else

  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float K=1*2000*2000;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;
  //printf("x_rank=%d,%d ,%d ,%d,%d,%d,%d\n",x_rank,stX,endX,stY,endY,stZ,endZ);
  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		K=vp_sub[n3d((i-5),(j-5),(k-5))]*vp_sub[n3d((i-5),(j-5),(k-5))];
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for(l=0;l< 5;l++){
		  temp1 += StaCoef[l] * (Vx[n3d_ext((i + l),j,k)] - Vx[n3d_ext((i - 1 - l),j,k)]);
		  temp2 += StaCoef[l] * (Vy[n3d_ext(i,(j+ l),k )] - Vy[n3d_ext(i,(j - 1 - l),k)]);
		  temp3 += StaCoef[l] * (Vz[n3d_ext(i,j,(k + l))] - Vz[n3d_ext(i,j,(k - 1 - l))]);
		}
		Txx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*Txx[n3d_ext(i,j,k)] + dt*K* (temp1) / dh);
		Tyy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*Tyy[n3d_ext(i,j,k)] + dt*K* (temp2) / dh);
		Tzz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*Tzz[n3d_ext(i,j,k)] + dt*K* (temp3) / dh);
	  }
	}
  }
#endif
}
/**************************
 * 反传传播相关函数*
 * ***********************/
//更新速度场
void update_Vel_back(){
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float rho=1.0f;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;

  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for(l=0;l < 5;l++){
		  temp1 += StaCoef[l] * 
			( Txx[n3d_ext((i+1+l),j,k)] - Txx[n3d_ext((i-l),j,k)]
			  +Tyy[n3d_ext((i+1+l),j,k)] - Tyy[n3d_ext((i-l),j,k)]
			  +Tzz[n3d_ext((i+1+l),j,k)] - Tzz[n3d_ext((i-l),j,k)]);

		  temp2 += StaCoef[l] * 
			(Txx[n3d_ext(i,(j+1+l),k)] - Txx[n3d_ext(i,(j-l),k)]
			 +Tyy[n3d_ext(i,(j+1+l),k)] - Tyy[n3d_ext(i,(j-l),k)]
			 +Tzz[n3d_ext(i,(j+1+l),k)] - Tzz[n3d_ext(i,(j-l),k)]);

		  temp3 += StaCoef[l] * 
			(Txx[n3d_ext(i,j,(k+1+l))] - Txx[n3d_ext(i,j,(k-l))]
			 +Tyy[n3d_ext(i,j,(k+1+l))] - Tyy[n3d_ext(i,j,(k-l))]
			 +Tzz[n3d_ext(i,j,(k+1+l))] - Tzz[n3d_ext(i,j,(k-l))]);
		}
		Vx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*Vx[n3d_ext(i,j,k)] - dt / rho* (temp1) / dh);
		Vy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*Vy[n3d_ext(i,j,k)] - dt / rho* (temp2) / dh);
		Vz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*Vz[n3d_ext(i,j,k)] - dt / rho * (temp3) / dh);
	  }
	}
  }
}
//更新应力场
void update_T_back(){
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float K=1*2000*2000;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;
  //printf("x_rank=%d,%d ,%d ,%d,%d,%d,%d\n",x_rank,stX,endX,stY,endY,stZ,endZ);
  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		K=vp_sub[n3d((i-5),(j-5),(k-5))]*vp_sub[n3d((i-5),(j-5),(k-5))];
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for (l=0;l< 5;l++){
		  temp1 += StaCoef[l] * (Vx[n3d_ext((i + l),j,k)] - Vx[n3d_ext((i - 1 - l),j,k)]);
		  temp2 += StaCoef[l] * (Vy[n3d_ext(i,(j+ l),k )] - Vy[n3d_ext(i,(j - 1 - l),k)]);
		  temp3 += StaCoef[l] * (Vz[n3d_ext(i,j,(k + l))] - Vz[n3d_ext(i,j,(k - 1 - l))]);
		}
		Txx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*Txx[n3d_ext(i,j,k)] - dt*K* (temp1) / dh);
		Tyy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*Tyy[n3d_ext(i,j,k)] - dt*K* (temp2) / dh);
		Tzz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*Tzz[n3d_ext(i,j,k)] - dt*K* (temp3) / dh);
	  }
	}
  }
}
//更新速度场
void update_Vel_extrap(){
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float rho=1.0f;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;

  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for(l=0;l < 5;l++){
		  temp1 += StaCoef[l] * 
			( extrapTxx[n3d_ext((i+1+l),j,k)] - extrapTxx[n3d_ext((i-l),j,k)]
			  +extrapTyy[n3d_ext((i+1+l),j,k)] - extrapTyy[n3d_ext((i-l),j,k)]
			  +extrapTzz[n3d_ext((i+1+l),j,k)] - extrapTzz[n3d_ext((i-l),j,k)]);

		  temp2 += StaCoef[l] * 
			(extrapTxx[n3d_ext(i,(j+1+l),k)] - extrapTxx[n3d_ext(i,(j-l),k)]
			 +extrapTyy[n3d_ext(i,(j+1+l),k)] - extrapTyy[n3d_ext(i,(j-l),k)]
			 +extrapTzz[n3d_ext(i,(j+1+l),k)] - extrapTzz[n3d_ext(i,(j-l),k)]);

		  temp3 += StaCoef[l] * 
			(extrapTxx[n3d_ext(i,j,(k+1+l))] - extrapTxx[n3d_ext(i,j,(k-l))]
			 +extrapTyy[n3d_ext(i,j,(k+1+l))] - extrapTyy[n3d_ext(i,j,(k-l))]
			 +extrapTzz[n3d_ext(i,j,(k+1+l))] - extrapTzz[n3d_ext(i,j,(k-l))]);
		}
		extrapVx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*extrapVx[n3d_ext(i,j,k)] - dt / rho* (temp1) / dh);
		extrapVy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*extrapVy[n3d_ext(i,j,k)] - dt / rho* (temp2) / dh);
		extrapVz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*extrapVz[n3d_ext(i,j,k)] - dt / rho * (temp3) / dh);
	  }
	}
  }
}
//更新应力场
void update_T_extrap(){
  //索引变量
  int i,j,k;
  //确定新的计算区域位置
  int stX=5,endX=SubNx-5;
  int stY=5,endY=SubNy-5;
  int stZ=5,endZ=SubNz-5;

  int X0_rank=5/DivNx;
  int Y0_rank=5/DivNy;
  int Z0_rank=5/DivNz;

  if(x_rank<X0_rank) stX=endX;
  if((dims[0]-(x_rank+1))<X0_rank) endX=stX;
  if(x_rank==X0_rank) stX=stX+(5-x_rank*DivNx);
  if((dims[0]-(x_rank+1))==X0_rank) endX=endX-(5-((dims[0]-(x_rank+1))*DivNx));

  if(y_rank<Y0_rank) stY=endY;
  if((dims[1]-(y_rank+1))<Y0_rank) endY=stY;
  if(y_rank==Y0_rank) stY=stY+(5-y_rank*DivNy);
  if((dims[1]-(y_rank+1))==Y0_rank) endY=endY-(5-((dims[1]-(y_rank+1))*DivNy));

  if(z_rank<Z0_rank) stZ=endZ;
  if((dims[2]-(z_rank+1))<Z0_rank) endZ=stZ;
  if(z_rank==Z0_rank) stZ=stZ+(5-z_rank*DivNz);
  if((dims[2]-(z_rank+1))==Z0_rank) endZ=endZ-(5-((dims[2]-(z_rank+1))*DivNz));

  int g_x,g_y,g_z;
  int l = 0;
  float ddx=0.0f,ddy=0.0f,ddz=0.0f;
  float K=1*2000*2000;
  float temp1 = 0.0f, temp2 = 0.0f,temp3=0.0f;
  //printf("x_rank=%d,%d ,%d ,%d,%d,%d,%d\n",x_rank,stX,endX,stY,endY,stZ,endZ);
  //Vx,Vy,Vz
  for(i=stX;i<endX;i++){
	g_x=x_rank*DivNx+i-5;
	ddx=pmlx_cal(g_x);
	for ( j = stY; j < endY; j++){
	  g_y=y_rank*DivNy+j-5;
	  ddy=pmly_cal(g_y);
	  for ( k = stZ; k < endZ; k++){
		K=vp_sub[n3d((i-5),(j-5),(k-5))]*vp_sub[n3d((i-5),(j-5),(k-5))];
		g_z=z_rank*DivNz+k-5;
		ddz=pmlz_cal(g_z);
		l=0;temp1 = 0.0, temp2 = 0.0, temp3 = 0.0;
		for (l=0;l< 5;l++){
		  temp1 += StaCoef[l] * (extrapVx[n3d_ext((i + l),j,k)] - extrapVx[n3d_ext((i - 1 - l),j,k)]);
		  temp2 += StaCoef[l] * (extrapVy[n3d_ext(i,(j+ l),k )] - extrapVy[n3d_ext(i,(j - 1 - l),k)]);
		  temp3 += StaCoef[l] * (extrapVz[n3d_ext(i,j,(k + l))] - extrapVz[n3d_ext(i,j,(k - 1 - l))]);
		}
		extrapTxx[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddx)*((1.0f - 0.5f*dt*ddx)*extrapTxx[n3d_ext(i,j,k)] - dt*K* (temp1) / dh);
		extrapTyy[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddy)*((1.0f - 0.5f*dt*ddy)*extrapTyy[n3d_ext(i,j,k)] - dt*K* (temp2) / dh);
		extrapTzz[n3d_ext(i,j,k)] = 1.0f / (1.0f + 0.5f*dt*ddz)*((1.0f - 0.5f*dt*ddz)*extrapTzz[n3d_ext(i,j,k)] - dt*K* (temp3) / dh);
	  }
	}
  }
}
